% posterior predictive simulation
clear

% catalog of data files
DFILE(1,:) = ['swuc_wpi_AR1_01'];
DFILE(2,:) = ['swuc_wpi_AR1_02'];
DFILE(3,:) = ['swuc_wpi_AR1_03'];
DFILE(4,:) = ['swuc_wpi_AR1_04'];
DFILE(5,:) = ['swuc_wpi_AR1_05'];
DFILE(6,:) = ['swuc_wpi_AR1_06'];
DFILE(7,:) = ['swuc_wpi_AR1_07'];
DFILE(8,:) = ['swuc_wpi_AR1_08'];
DFILE(9,:) = ['swuc_wpi_AR1_09'];
DFILE(10,:) = ['swuc_wpi_AR1_10'];
DFILE(11,:) = ['swuc_wpi_AR1_11'];
DFILE(12,:) = ['swuc_wpi_AR1_12'];
DFILE(13,:) = ['swuc_wpi_AR1_13'];
DFILE(14,:) = ['swuc_wpi_AR1_14'];
DFILE(15,:) = ['swuc_wpi_AR1_15'];
DFILE(16,:) = ['swuc_wpi_AR1_16'];
DFILE(17,:) = ['swuc_wpi_AR1_17'];
DFILE(18,:) = ['swuc_wpi_AR1_18'];
DFILE(19,:) = ['swuc_wpi_AR1_19'];
DFILE(20,:) = ['swuc_wpi_AR1_20'];
DFILE(21,:) = ['swuc_wpi_AR1_21'];
DFILE(22,:) = ['swuc_wpi_AR1_22'];
DFILE(23,:) = ['swuc_wpi_AR1_23'];
DFILE(24,:) = ['swuc_wpi_AR1_24'];
DFILE(25,:) = ['swuc_wpi_AR1_25'];
DFILE(26,:) = ['swuc_wpi_AR1_26'];
DFILE(27,:) = ['swuc_wpi_AR1_27'];
DFILE(28,:) = ['swuc_wpi_AR1_28'];
DFILE(29,:) = ['swuc_wpi_AR1_29'];
DFILE(30,:) = ['swuc_wpi_AR1_30'];
DFILE(31,:) = ['swuc_wpi_AR1_31'];
DFILE(32,:) = ['swuc_wpi_AR1_32'];
DFILE(33,:) = ['swuc_wpi_AR1_33'];
DFILE(34,:) = ['swuc_wpi_AR1_34'];
DFILE(35,:) = ['swuc_wpi_AR1_35'];
DFILE(36,:) = ['swuc_wpi_AR1_36'];
DFILE(37,:) = ['swuc_wpi_AR1_37'];
DFILE(38,:) = ['swuc_wpi_AR1_38'];
DFILE(39,:) = ['swuc_wpi_AR1_39'];
DFILE(40,:) = ['swuc_wpi_AR1_40'];

NF = size(DFILE,1);
NB = 21;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read data

% warren-pearson data 1749-1890 (1782-1784, 1788, 1792 are missing) 1914 = 100 
% Historical Statistics of the United States
% read data for 1797-1890
wp = xlsread('c:\wp inflation volatility\US\wholesale\Cc113-124','B45:B138'); 
inf_wp = diff(log(wp));
T_wp = size(wp,1);
year_wp = 1797 + [0:1:T_wp-1];
% figure
% plot(year_wp(2:T_wp),inf_wp,'-b','Linewidth',2); hold on;
% axis([1795 2015 -inf inf])

% hanes data 1860-1990 (1942-46 are missing), 
% Historical Statistics of the United States
% For the period 1860-1941, 1914 = 100
hanes1 = xlsread('c:\wp inflation volatility\US\wholesale\Cc125-137','B3:B84'); 
inf_h1 = diff(log(hanes1));
T_hanes1 = size(hanes1,1);
year_hanes1 = 1860 + [0:1:T_hanes1-1];
% plot(year_hanes1(2:T_hanes1),inf_h1,'-g','Linewidth',2); hold on;

% For the period 1947-1990, 1990 = 100
hanes2 = xlsread('c:\wp inflation volatility\US\wholesale\Cc125-137','B86:B129'); 
inf_h2 = diff(log(hanes2));
T_hanes2 = size(hanes2,1);
year_hanes2 = 1947 + [0:1:T_hanes2-1];
% plot(year_hanes2(2:T_hanes2),inf_h2,'-g','Linewidth',2); hold on;

% bls data 1890-1997 (1982 = 100)
% Historical Statistics of the United States
bls_hs = xlsread('c:\wp inflation volatility\US\wholesale\Cc66-83','B2:B109'); 
inf_bls_hs = diff(log(bls_hs));
T_bls_hs = size(bls_hs,1);
year_bls_hs = 1890 + [0:1:T_bls_hs-1];
%plot(year_bls_hs(2:T_bls_hs),inf_bls_hs,'-r','Linewidth',2); hold on;

% bls data from FRED (1982 = 100), monthly, NSA, 1/1913-6/2013
bls_m_fred = xlsread('c:\wp inflation volatility\US\wholesale\PPI_ACO_FRED_July_16_2013','B13:B1218'); % monthly
Tm = size(bls_m_fred,1);
bls_fred = bls_m_fred(6:12:Tm); % annual, point sampled in June (maximizes correlation with annual BLS HSUS data)
inf_bls_fred = diff(log(bls_fred));
T_bls_fred = size(bls_fred,1);
year_bls_fred = 1913 + [0:1:T_bls_fred-1];
% plot(year_bls_fred(2:T_bls_fred),inf_bls_fred,'-m','Linewidth',2); hold off;
% legend('Warren-Pearson','Hanes 1891-1941','Hanes 1948-1990','BLS HSUS','BLS FRED')
% set(gca,'Fontsize',16)

% noisy inflation data Warren-Pearson 1798-1860,Hanes 1861-1941, BLS HSUS
% 1942-1947, Hanes 1948-1990
yn = [inf_wp(1:63); inf_h1; inf_bls_hs(52:57); inf_h2]; % noisy inflation data 
% clean inflation data
yc = inf_bls_fred(35:T_bls_fred-1); % BLS FRED 1948-2013
year = [1798:1:2013]';
T = size(year,1);
% figure
% plot(year(1:size(yn,1)),yn,'-g','Linewidth',2'); hold on;
% plot(year(151:T),yc,'-b','Linewidth',2'); hold off;
% legend('Noisy','Clean')
% set(gca,'Fontsize',16)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% partition sample
Y0 = yn(1:52); % training sample 1798-1849

year = [1850:1:2013]'; % years to be used for estimation
T1 = 98; % 1947
T2 = 141; %1990
T = size(year,1);

YS1 = yn(53:size(yn,1))'; % 1850-1990 to be used for estimation
YS2 = [zeros(1,T1) yc']; % 1948-2013 padded w/ zeros to make dating easier

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% standard deviation of innovations to log volatilities
% load WA
load(DFILE(NB,:),'WD')
[P,N,N] = size(WD);
WA = zeros((NF-NB+1)*P,N,N);
WA(1:P,:,:) = WD;
clear WD 
for i = NB+1:NF,
  j = i - NB + 1; 
  load(DFILE(i,:),'WD')
  WA((j-1)*P+1:j*P,:,:) = WD;
  clear WD;
end
NMC = size(WA,1);
W11 = squeeze(WA(:,1,1));
W22 = squeeze(WA(:,2,2));
W12 = squeeze(WA(:,1,2));
clear WA

% prior for W (inverse wishart w/ scale matrix W0 and degrees of freedom v0) (innovation variance for log volatilities)
v0 = 10;
sv0 = 0.2236*sqrt((v0+1)/v0); % stock and watson's calibrated value adjusted for time aggregation
TW0 = v0*(sv0^2)*eye(2);
PS = real(sqrtm(TW0\eye(2,2))); 
Wprior = zeros(100000,2,2);
for ii = 1:100000,
    u = randn(2,v0);
    Wprior(ii,:,:) = (PS*u*u'*PS')\eye(2,2);
end
W11p = squeeze(Wprior(:,1,1));
W22p = squeeze(Wprior(:,2,2));
W12p = squeeze(Wprior(:,1,2));
clear Wprior 

grid = 0.1:.01:0.75;
[Np1,Xp1] = hist(W11p.^.5,grid);
[Nr,Xr] = hist(W11.^.5,grid);
[Nq,Xq] = hist(W22.^.5,grid);

figure
subplot(2,2,1)
plot(Xp1,Np1/sum(Np1'),'--k','Linewidth',2); hold on;
plot(Xr,Nr/sum(Nr'),'-k','Linewidth',2); hold on;
plot(Xq,Nq/sum(Nq'),':k','Linewidth',2); hold off;
xlabel('\sigma_r,\sigma_q','Fontsize',18)
legend('Prior','Posterior \sigma_r','Posterior \sigma_q')
axis([0 0.7 0 inf])
set(gca,'Fontsize',16)

corr_prior = W12p./((W11p.*W22p).^.5);
corr_post = W12./((W11.*W22).^.5);
clear W11 W22 W12 W11p W22p W12p 

grid = -1:.05:1;
[Nprq,Xprq] = hist(corr_prior,grid);
[Nrq,Xrq] = hist(corr_post,grid);
subplot(2,2,2)
plot(Xprq,Nprq/sum(Nprq'),'--k','Linewidth',2); hold on;
plot(Xrq,Nrq/sum(Nrq'),'-k','Linewidth',2); hold off;
xlabel('Correlation','Fontsize',18)
legend('Prior','Posterior')
axis([-1 1 0 inf])
set(gca,'Fontsize',16)

clear Xp1 Np1 Xr Nr Xq Nq Xprq Nprq Xnq Nrq 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% standard deviation of pre-1948 measurement errors
% load MA
load(DFILE(NB,:),'MD')
[P,N] = size(MD);
MA = zeros(N,(NF-NB+1)*P);
MA(:,1:P) = MD';
clear MD 
for i = NB+1:NF,
  j = i - NB + 1; 
  load(DFILE(i,:),'MD')
  MA(:,(j-1)*P+1:j*P) = MD';
  clear MD;
end

% prior for measurement-error variance \sigma_m (prior is same for all periods)
vm0 = 2;
R0 = .0833^2;
sm0 = 0.5*sqrt(R0)*sqrt((vm0+1)/vm0);
dm0 = vm0*(sm0^2);
svmc2 = zeros(100000,1);
for ii = 1:100000,
    v = ig2(vm0,dm0,0);
    svmc2(ii,1) = v^.5;
end

grid = [.015:.005:.26];
[Np2,Xp2] = hist(svmc2(:,1),grid);
[Nm1,Xm1] = hist(MA(1,:),grid);

figure
subplot(2,2,2)
plot(Xp2,Np2/sum(Np2'),'--k','Linewidth',2); hold on;
plot(Xm1,Nm1/sum(Nm1'),'-k','Linewidth',2); hold off;
xlabel('\sigma_{m}','Fontsize',18)
%title('Prior and Posteriors for the Standard Deviation of Measurement Innovations','Fontsize',16)
legend('Prior \sigma_m','Posterior \sigma_m')
axis([0 0.15 0 inf])
set(gca,'Fontsize',16)

clear MA svmc2 Xp2 Np2 Xm1 Nm1 grid

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prior for measurement-error ar parameter \rho 
rho0 = 0; % prior mean on ar coefficients
Vrho0 = 0.45^2; % prior variance on ar coefficients
invVrho0 = inv(Vrho0);

% autoregressive parameters
% load AA
load(DFILE(NB,:),'AD')
[P,N] = size(AD);
AA = zeros(N,(NF-NB+1)*P);
AA(:,1:P) = AD';
clear AD 
for i = NB+1:NF,
  j = i - NB + 1; 
  load(DFILE(i,:),'AD')
  AA(:,(j-1)*P+1:j*P) = AD';
  clear AD;
end
NMC = size(AA,2);
smc_rho = ones(NMC,1)*rho0' + (sqrtm(Vrho0)*randn(1,NMC))'; %prior for rho
[Np1,Xp1] = hist(smc_rho(:,1),50); % prior
[Na1,Xa1] = hist(AA(1,:),Xp1);
subplot(2,2,1)
plot(Xp1,Np1/sum(Np1),'--k','Linewidth',2); hold on;
plot(Xa1,Na1/sum(Na1),'-k','Linewidth',2); hold off;
legend('Prior \rho_m','Posterior \rho_m')
xlabel('\rho_{m}','Fontsize',18)
axis([-1 1 0 inf])
set(gca,'Fontsize',16)

'Probability that \rho exceeds 0'
DD = zeros(size(AA));
DD(AA>0) = 1.0;
mean(DD)

'Probability that \rho is explosive'
DD = zeros(size(AA));
DD(abs(AA)>1) = 1.0;
mean(DD)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  compute RA and QA
load(DFILE(NB,:),'vD')
[N,T,P] = size(vD);
RA = zeros(T,(NF-NB+1)*P);
QA = zeros(T,(NF-NB+1)*P);
RA(:,1:P) = vD(1,:,:);
QA(:,1:P) = vD(2,:,:);
clear vD 
for i = NB+1:NF,
  j = i - NB + 1; 
  load(DFILE(i,:),'vD')
  RA(:,(j-1)*P+1:j*P) = vD(1,:,:);
  QA(:,(j-1)*P+1:j*P) = vD(2,:,:);
  clear vD;
end
SRA = sort(RA(2:T,:),2);
NMC = size(SRA,2);
CRA = SRA(:,NMC*[.25 .5 .75]); % median and interquartile range for RA
clear RA SRA
SQA = sort(QA(2:T,:),2);
CQA = SQA(:,NMC*[.25 .5 .75]); % median and interquartile range for QA
clear QA SQA

% compute SA
load(DFILE(NB,:),'SD')
[P,K,T] = size(SD);
SA1 = zeros((NF-NB+1)*P,T); % \pi
SA2 = zeros((NF-NB+1)*P,T); % \mu
SA1(1:P,:,:) = squeeze(SD(:,1,:));
SA2(1:P,:,:) = squeeze(SD(:,2,:));
clear SD 
for i = NB+1:NF,
  j = i - NB + 1; 
  load(DFILE(i,:),'SD')
  SA1((j-1)*P+1:j*P,:,:) = squeeze(SD(:,1,:));
  SA2((j-1)*P+1:j*P,:,:) = squeeze(SD(:,2,:));
  clear SD;
end
SSA2 = sort(SA2,1)';
CSA2 = SSA2(:,NMC*[.25 .5 .75]); % median and interquatile range 
clear SSA2

SAT = SA1-SA2; % transient component
clear SA1 SA2

SSAT = sort(SAT,1)';
CSAT = SSAT(:,NMC*[.25 .5 .75]); % median and interquatile range for RA
clear SAT SSAT

figure
subplot(2,2,1)
plot(year,CSA2(:,2),'-k','Linewidth',2); hold on;
plot(year,CSA2(:,1),':k','Linewidth',2); hold on;
plot(year,CSA2(:,3),':k','Linewidth',2); hold off;
title('\mu_{t}','Fontsize',16)
legend('Median','Interquartile range','Location','Southeast')
set(gca,'Fontsize',14)
axis([1850 2015 -inf inf])
grid

subplot(2,2,2)
plot(year(1:5:T),CSA2(1:5:T,2),'--k','Linewidth',2); hold on;
plot(year,CSAT(:,2),'-k','Linewidth',2); hold on;
plot(year,CSAT(:,1),':k','Linewidth',2); hold on;
plot(year,CSAT(:,3),':k','Linewidth',2); hold off;
title('\pi_{t} - \mu_{t}','Fontsize',16)
set(gca,'Fontsize',14)
axis([1850 2015 -inf inf])
grid

subplot(2,2,3)
plot(year,CQA(:,2).^.5,'-k','Linewidth',2); hold on;
plot(year,CQA(:,1).^.5,':k','Linewidth',2); hold on;
plot(year,CQA(:,3).^.5,':k','Linewidth',2); hold off;
title('q_{t}^{1/2}','Fontsize',16)
set(gca,'Fontsize',14)
axis([1850 2015 0 inf])
grid

subplot(2,2,4)
plot(year,CRA(:,2).^.5,'-k','Linewidth',2); hold on;
plot(year,CRA(:,1).^.5,':k','Linewidth',2); hold on;
plot(year,CRA(:,3).^.5,':k','Linewidth',2); hold off;
title('r_{t}^{1/2}','Fontsize',16)
set(gca,'Fontsize',14)
axis([1850 2015 0 inf])
grid

